
# Import Libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import plotly
import plotly.offline as pyoff
import plotly.graph_objs as go
import plotly.express as px
import chart_studio
import chart_studio.plotly as py
import calmap
import datetime
import tensorflow as tf
import os
import random
import re
import plotly.offline as pyoff
import plotly.graph_objs as go

from datetime import date
from plotly.subplots import make_subplots
from itertools import cycle, product
from statsmodels.tsa.seasonal import STL
from scipy.stats import boxcox
from pmdarima.arima import ADFTest, auto_arima
from pmdarima.utils import diff_inv
from statsmodels.tsa.stattools import adfuller
from sklearn.model_selection import TimeSeriesSplit
from tensorflow.keras.layers import LSTM, Dense, BatchNormalization
from tensorflow.keras import Sequential
from tensorflow.keras.backend import clear_session
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator
from tensorflow.keras.initializers import *
from tensorflow.keras import optimizers
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.linear_model import LinearRegression
from scipy.special import boxcox1p, inv_boxcox1p
import matplotlib.patches as mpatches
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.model_selection import GridSearchCV
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.optimizers import Adam
from IPython.core.display import HTML, display


# Versões dos pacotes usados neste jupyter notebook
%reload_ext watermark
%watermark -a "Herikc Brecher" --iversions

Author: Herikc Brecher

pandas      : 1.2.2
tensorflow  : 2.5.0
calmap      : 0.0.9
re          : 2.2.1
numpy       : 1.19.5
plotly      : 4.14.3
datetime    : 4.3
chart_studio: 1.1.0
matplotlib  : 3.4.2
seaborn     : 0.11.1


# Variaveis globais
SEED = 84796315
FEATURES = 7
EPOCHS = 100
BATCH_SIZE = 1000
EXECUTE_GRID_SEARCH = True


# Configurando seeds
os.environ['PYTHONHASHSEED'] = str(SEED)
tf.random.set_seed(SEED)
np.random.seed(SEED)
random.seed(SEED)


# Exibindo toda tela
display(HTML('<style>.container { width:100% !important; }</style>'))
pd.options.plotting.backend = 'matplotlib'


# Import dataset
dtOrders = pd.read_csv('../data/olist_orders_dataset.csv', encoding = 'utf8')


# Colunas do tipo data
dateColumns = ['order_purchase_timestamp', 'order_approved_at', 'order_delivered_carrier_date',\
               'order_delivered_customer_date', 'order_estimated_delivery_date']

# Dataset de analise temporal
dtOrdersAdjusted = dtOrders.copy()


# Convertendo columas de data para date
for col in dateColumns:
    dtOrdersAdjusted[col] = pd.to_datetime(dtOrdersAdjusted[col], format = '%Y-%m-%d %H:%M:%S')


# Dropando valores NA
dtOrdersAdjusted = dtOrdersAdjusted.dropna()


dtOrdersAdjusted.dtypes

order_id                                 object
customer_id                              object
order_status                             object
order_purchase_timestamp         datetime64[ns]
order_approved_at                datetime64[ns]
order_delivered_carrier_date     datetime64[ns]
order_delivered_customer_date    datetime64[ns]
order_estimated_delivery_date    datetime64[ns]
dtype: object


dtHistory = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date

start = dtHistory.min()
end = dtHistory.max()

idx = pd.date_range(start, end, normalize = True)

seriesHistory = dtHistory.value_counts(sort = False).sort_index().reindex(idx, fill_value = 0)

dtHistory = pd.DataFrame(seriesHistory).reset_index()


dtHistory


# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['index'],
                        y = dtHistory['order_purchase_timestamp'])]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {'title': 'Vendas'},
                        title = 'Vendas por dia')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

pyoff.iplot(fig)


# Remove outliers
seriesHistory = seriesHistory[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
pred_range = pd.date_range(datetime.date(2018, 8, 17), datetime.date(2018, 10, 17))
dtHistory = pd.DataFrame(seriesHistory).reset_index()


# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['index'],
                        y = dtHistory['order_purchase_timestamp'])]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {'title': 'Vendas'},
                        title = 'Vendas por dia')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

pyoff.iplot(fig)


#Plot histórico de vendas por dia, mês e ano
fig, caxs = calmap.calendarplot(seriesHistory, daylabels = 'MTWTFSS', fillcolor = 'grey',cmap = 'YlGn', fig_kws = dict(figsize = (18, 9)))
fig.suptitle('Histórico de Vendas', fontsize = 22)

fig.subplots_adjust(right = 0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.67])
fig.colorbar(caxs[0].get_children()[1], cax = cbar_ax)

plt.show()


# Criar grafico na estrutura STL 4 layers
def add_stl_plot(fig, res, legend):
    axs = fig.get_axes()
    
    # Nome de cada um dos subplots
    comps = ['trend', 'seasonal', 'resid']
    for ax, comp in zip(axs[1:], comps):
        series = getattr(res, comp)
        if comp == 'resid':
            ax.plot(series, marker = 'o', linestyle = 'none')
        else:
            ax.plot(series)
            ax.legend(legend, frameon = False)


# Gerar STL
stl = STL(seriesHistory)
stl_res = stl.fit()
fig = stl_res.plot()
fig.set_size_inches((20, 12))
plt.show()


# Gerar STL não robusto e concatenar ao robusto
stl = STL(seriesHistory, robust = True)
res_robust = stl.fit()
fig = res_robust.plot()
fig.set_size_inches((20, 12))
res_non_robust = STL(seriesHistory, robust = False).fit()
add_stl_plot(fig, res_non_robust, ['Robusto', 'Não Robusto'])


# Gerando STL para separar cada um dos componentes
stl = STL(seriesHistory)
res = stl.fit()

# Separando seriesDeseasonal
seriesDeseasonal = res.observed - res.seasonal

# Separando boxcox
seriesBoxCox, lmbda = boxcox(seriesHistory+1)
seriesBoxCox = pd.Series(seriesBoxCox, index = seriesHistory.index)

# Separando stationary
seriesDiffHistory = seriesHistory.diff(7).dropna()
xi = seriesHistory.iloc[:7]

# Plotando todos elementos de forma separada
fig, axs = plt.subplots(nrows = 4, sharex = True, figsize = (14, 8))

seriesHistory.plot(ax = axs[0])
axs[0].set_ylabel('Original')

seriesDeseasonal.plot(ax = axs[1])
axs[1].set_ylabel('Deseasonal')

seriesBoxCox.plot(ax = axs[2])
axs[2].set_ylabel('BoxCox')

seriesDiffHistory.plot(ax = axs[3])
axs[3].set_ylabel('Stationary')

fig.align_ylabels()
fig.suptitle('Transformações')

plt.tight_layout()
plt.show()


# Teste ADF automatizado
def adf_test(series):
    result = adfuller(series, autolag = 'AIC')

    print(f'ADF: {result[0]}')
    print(f'steps: {result[1]}')
    print(f'p-value: {result[1]}')

    for key, value in result[4].items():
        print('\nValor Critico:')
        print(f'   {key}: {value}')


adf = ADFTest(alpha = 0.05)
adf.should_diff(seriesHistory)

(0.01, False)


# Testando estacionaridade na série original
adf_test(seriesHistory)

ADF: -2.6163082378564684
steps: 0.08967331201291256
p-value: 0.08967331201291256

Valor Critico:
   1%: -3.441694608475642

Valor Critico:
   5%: -2.866544718556839

Valor Critico:
   10%: -2.5694353738653684


adf.should_diff(seriesDiffHistory)

(0.01, False)


adf_test(seriesDiffHistory)

ADF: -6.801606152090551
steps: 2.229745766641977e-09
p-value: 2.229745766641977e-09

Valor Critico:
   1%: -3.441834071558759

Valor Critico:
   5%: -2.8666061267054626

Valor Critico:
   10%: -2.569468095872659


# Controle de resultados de toda fase de modelagem
result = pd.DataFrame(columns = ['Algorithm', 'MSE', 'RMSE', 'MAE', 'Mean_Real_Value', 'Mean_Predict_Value'])


split_range = TimeSeriesSplit(n_splits = 8, max_train_size = pred_range.shape[0], test_size = pred_range.shape[0])


# Adiciona o registro ao dataset
def record(result, algorithm, mse = -1, rmse = -1, mae = -1, mrv = -1, mpv = -1, show = True):
    new = pd.DataFrame(dict(Algorithm = algorithm, MSE = mse, RMSE = rmse, MAE = mae, Mean_Real_Value = mrv,\
                            Mean_Predict_Value = mpv), index = [0])
    result = pd.concat([result, new], ignore_index = True)
    
    if show:
        display(result)
    
    return result


# Plot no formato de 4 layers, seguindo o STL para cada um dos modelos
def plot(index, pred, mse, title, fig = None, ax = None, ylabel = ''):
    global seriesHistory
    
    empty_fig = fig is None
    
    if empty_fig:
        fig, ax = plt.subplots(figsize = (13, 6))
    else: 
        ax.set_ylabel(ylabel)
                
    ax.set_title(title)    
    patch_ = mpatches.Patch(color = 'white', label = f'MSE: {np.mean(mse):.1e}')
    L1 = ax.legend(handles = [patch_], loc = 'upper left', fancybox = True, framealpha = 0.7,  handlelength = 0)
    ax.add_artist(L1)
    
    sns.lineplot(x = seriesHistory.index, y = seriesHistory, label = 'Real', ax = ax)
    sns.lineplot(x = index, y = pred, label = 'Previsto', ax = ax)
    ax.axvline(x = index[0], color = 'red')
    ax.legend(loc = 'upper right')
    
    if empty_fig:
        plt.show()
    else:
        return fig


# Calculo para previsão e teste quando utilizado a série Original
def calcPredTestOriginal(train, pred, test):
    return pred, test, 0


# Calculo para previsão e teste quando utilizado a série seriesDeseasonal
def calcPredTestseriesDeseasonal(train, pred, test):
    # Removendo a sazonalidade da série e convertendo para o shape correto
    last_seasonal = res.seasonal.reindex_like(train).tail(stl.period)
    pred = pred + np.fromiter(cycle(last_seasonal), count = pred.shape[0], dtype = float)
    test = test + res.seasonal.reindex_like(test)
    
    return pred, test, 1


# Calculo para previsão e teste quando utilizado a série BoxCox
def calcPredTestBoxCox(train, pred, test):
    # Reverdendo a normalização do boxcox
    pred = inv_boxcox1p(pred, lmbda)
    test = inv_boxcox1p(test, lmbda)
    
    return pred, test, 2


# Calculo para previsão e teste quando utilizado a série Stationary
def calcPredTestStationary(train, pred, test):
    # Calculando a diferença da sazonalidade
    xi = seriesHistory.reindex_like(train).tail(FEATURES)
    
    totalLen = len(pred) + len(xi) 
    ix = pd.date_range(xi.index[0], periods = totalLen)  
    inv = diff_inv(pred, FEATURES, xi = xi) + np.fromiter(cycle(xi), count = totalLen, dtype = float)  
    inv = pd.Series(inv, index = ix, name = 'order_purchase_timestamp')
    pred = inv.iloc[FEATURES:]
    
    totalLen = len(test) + len(xi) 
    ix = pd.date_range(xi.index[0], periods = totalLen)  
    inv = diff_inv(test, FEATURES, xi = xi) + np.fromiter(cycle(xi), count = totalLen, dtype = float)  
    inv = pd.Series(inv, index = ix, name = 'order_purchase_timestamp')
    test = inv.iloc[FEATURES:]
    
    return pred, test, 3


# Report para Time Series Regressor, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportTSR(data, modelName, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - Time Series Regression'
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
    
        gen = TimeseriesGenerator(train, train, FEATURES, batch_size = BATCH_SIZE)

        X_train = gen[0][0]
        y_train = gen[0][1]

        lr = LinearRegression()
        lr.fit(X_train, y_train)
        X_pred = y_train[-FEATURES:].reshape(1,-1)
        pred = np.empty(test.shape[0])

        for i in range(len(pred)):
            forecast = lr.predict(X_pred)
            X_pred = np.delete(X_pred, 0, 1)
            X_pred = np.concatenate((X_pred, forecast.reshape(-1, 1)), 1)    
            pred[i] = forecast
        
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.tight_layout()
plt.close()


reportTSR(seriesHistory.copy(), 'Original', calcPredTestOriginal)


reportTSR(seriesDeseasonal.copy(), 'Deseasonal', calcPredTestseriesDeseasonal)


reportTSR(seriesBoxCox.copy(), 'BoxCox', calcPredTestBoxCox)


reportTSR(seriesDiffHistory.copy(), 'Stationary', calcPredTestStationary)


result


# Função utilizada para o hypertuning de alpha, beta e gamma do Exponential Smoothing
def GSES(data, modelName, alpha, beta, gamma, calcFunction):    
    mse = []
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        try:
            with catch_warnings():
                filterwarnings('ignore')
                ES = (
                    ExponentialSmoothing(train, trend = 'add', seasonal = 'add', seasonal_periods = FEATURES)
                    .fit(smoothing_level = alpha, smoothing_trend = beta, smoothing_seasonal = gamma, method = 'ls')
                )

                pred = ES.forecast(test.shape[0])

                pred, test, _ = calcFunction(train, pred, test)

                mse.append(mean_squared_error(pred, test, squared = True))
        
        except:
            mse.append(-1)
    
    return np.mean(mse)


# Função utilizada para o hypertuning de demais parâmetros do Exponential Smoothing
def GSESOPT(data, modelName, trend, season, periods, bias, method, calcFunction):
    mse = []
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        try:
            with catch_warnings():
                filterwarnings('ignore')
                ES = (
                    ExponentialSmoothing(train, trend = trend, seasonal = season, seasonal_periods = periods)
                    .fit(remove_bias = bias, method = method, optimized = True)
                )

                pred = ES.forecast(test.shape[0])

                pred, test, _ = calcFunction(train, pred, test)

                mse.append(mean_squared_error(pred, test, squared = True))       
        except:
            mse.append(-1)
    
    return np.mean(mse)


# Report para Exponential Smoothing, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportES(data, modelName, model_kwargs, fit_kwargs, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - Exponential Smoothing'
    indexPlot = 0
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        ES = ExponentialSmoothing       
        stlf = STLForecast(train, ES, model_kwargs = model_kwargs).fit(fit_kwargs = fit_kwargs)
        pred = stlf.forecast(test.shape[0])
             
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Função para gerar tabela de hypertuning ampla
def exp_smoothing_configs(seasonal = [None]):
    models = list()
    # Lista de argumentos
    t_params = ['add', 'mul']
    s_params = ['add', 'mul']
    p_params = seasonal
    r_params = [True, False]
    method_params = ['L-BFGS-B' , 'TNC', 'SLSQP', 'Powell', 'trust-constr', 'bh', 'ls']
    
    # Gerando lista de argumentos
    for t in t_params:
        for s in s_params:
            for p in p_params:
                for r in r_params:
                    for m in method_params:
                        cfg = [t, s, p, r, m]
                        models.append(cfg)
    return models


# Gerando tabela de hypertunning
alphas = betas = gammas = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas), columns = ['alpha', 'beta', 'gamma'])
hyperparam.head()


%%time

# Treinamento do modelo 
if EXECUTE_GRID_SEARCH:  
    hyperparam['mse'] = hyperparam.apply(lambda x: GSES(seriesHistory.copy(), 'Original',\
                                                x.alpha, x.beta, x.gamma, calcPredTestOriginal), axis = 1)

Wall time: 59.2 s


# Verificando o menor mse
if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


# Criando lista de argumentos ampla
params_ = exp_smoothing_configs([FEATURES])


hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])


len(hyperparam_)

56


hyperparam_.head()


%%time

# Se True irá treinar com a nova lista mais ampla (pode demorar)
if EXECUTE_GRID_SEARCH:
    hyperparam_['mse'] = hyperparam_.apply(lambda x: GSESOPT(seriesHistory.copy(), 'Original',\
                                                     x.trend, x.season, x.periods, x.bias, x.method, calcPredTestOriginal),\
                                   axis = 1)

Wall time: 3min 21s


if EXECUTE_GRID_SEARCH:
    display(hyperparam_.query('mse == mse.min() and mse != -1'))


# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13, 6))
figs.align_ylabels()
figs.tight_layout()
plt.close()


# Dicionar de configurações de acordo com o hypertuning
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES)
fit_kwargs = dict(smoothing_level = 0.2, smoothing_trend = 0, smoothing_seasonal = 0.8, method = 'ls')


reportES(seriesHistory.copy(), 'Original', model_kwargs, fit_kwargs, calcPredTestOriginal)


alphas = betas = gammas = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas), columns = ['alpha', 'beta', 'gamma'])
hyperparam.head()


%%time

if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSES(seriesDeseasonal.copy(), 'seriesDeseasonal',\
                                                x.alpha, x.beta, x.gamma, calcPredTestseriesDeseasonal), axis = 1)

Wall time: 1min 2s


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


params_ = exp_smoothing_configs([FEATURES])


hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])


%%time
if EXECUTE_GRID_SEARCH:
    hyperparam_['mse'] = hyperparam_.apply(lambda x: GSESOPT(seriesDeseasonal.copy(), 'seriesDeseasonal',\
                                                     x.trend, x.season, x.periods, x.bias, x.method, calcPredTestseriesDeseasonal),\
                                   axis = 1)

Wall time: 4min 21s


if EXECUTE_GRID_SEARCH:
    display(hyperparam_.query('mse == mse.min() and mse != -1'))


model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES)
fit_kwargs = dict(smoothing_level = 0.2, smoothing_trend = 0, smoothing_seasonal = 0.8, method = 'ls')


reportES(seriesDeseasonal.copy(), 'Deseasonal', model_kwargs, fit_kwargs, calcPredTestseriesDeseasonal)


alphas = betas = gammas = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas), columns = ['alpha', 'beta', 'gamma'])
hyperparam.head()


%%time

if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSES(seriesBoxCox.copy(), 'BoxCox',\
                                                x.alpha, x.beta, x.gamma, calcPredTestBoxCox), axis = 1)

Wall time: 1min


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


params_ = exp_smoothing_configs([FEATURES])


hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])


%%time
if EXECUTE_GRID_SEARCH:
    hyperparam_['mse'] = hyperparam_.apply(lambda x: GSESOPT(seriesBoxCox.copy(), 'BoxCox',\
                                                  x.trend, x.season, x.periods, x.bias, x.method, calcPredTestBoxCox),\
                                   axis = 1)

Wall time: 3min 43s


if EXECUTE_GRID_SEARCH:
    display(hyperparam_.query('mse == mse.min() and mse != -1'))


model_kwargs = dict(trend = 'mul', seasonal = 'add', seasonal_periods = FEATURES)
fit_kwargs = dict(remove_bias = True, smoothing_level = 0.7, smoothing_trend = 0.2, smoothing_seasonal = 0.8, method = 'ls')


reportES(seriesBoxCox.copy(), 'BoxCox', model_kwargs, fit_kwargs, calcPredTestBoxCox)


alphas = betas = gammas = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas), columns = ['alpha', 'beta', 'gamma'])
hyperparam.head()


%%time
if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSES(seriesDiffHistory.copy(), 'Stationary',\
                                                x.alpha, x.beta, x.gamma, calcPredTestStationary), axis = 1)

Wall time: 1min 7s


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


params_ = exp_smoothing_configs([FEATURES])


hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])


%%time
if EXECUTE_GRID_SEARCH:
    hyperparam_['mse'] = hyperparam_.apply(lambda x: GSESOPT(seriesDiffHistory.copy(), 'Stationary',\
                                                  x.trend, x.season, x.periods, x.bias, x.method, calcPredTestStationary),\
                                   axis = 1)

Wall time: 1min 37s


if EXECUTE_GRID_SEARCH:
    display(hyperparam_.query('mse == mse.min() and mse != -1'))


model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES)
fit_kwargs = dict(smoothing_level = 0, smoothing_trend = 0.3, smoothing_seasonal = 0, method = 'ls')


reportES(seriesDiffHistory.copy(), 'Stationary', model_kwargs, fit_kwargs, calcPredTestStationary)


# Report do algoritmo arima, também é adicionado a base de resultados e realizado o plot de acompanhamento
def reportArima(arimaModel, modelName, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - '  + arimaModel.__str__().strip()
    indexPlot = 0
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        arimaModel.fit(train)
        pred = arimaModel.predict(test.shape[0])
             
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()


# Correlação entre os periodos com ARIMA

lags = 90

with catch_warnings():
    filterwarnings('ignore')
    fig, ax = plt.subplots(2, figsize = (12, 6), sharex = True)
    plot_acf(seriesHistory.diff().dropna(), ax = ax[0], lags = lags, missing = 'drop')
    plot_pacf(seriesHistory.diff().dropna(), ax = ax[1], lags = lags)
    plt.show()


%%time

# Utilizando o auto arima para descobrir os argumentos ideias baseados no conjunto de dado informado
data = seriesHistory.copy()
arimaModel = auto_arima(seriesHistory.copy(), m = FEATURES, seasonal = True)
arimaModel

Wall time: 22.2 s

ARIMA(order=(1, 1, 2), scoring_args={}, seasonal_order=(0, 0, 2, 7),
      suppress_warnings=True)


reportArima(arimaModel, 'Original', calcPredTestOriginal)


%%time
data = seriesDeseasonal.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = False)
arimaModel

C:\Users\herik\AppData\Roaming\Python\Python38\site-packages\pmdarima\arima\_validation.py:62: UserWarning:

m (7) set for non-seasonal fit. Setting to 0

Wall time: 2.6 s

ARIMA(order=(2, 1, 1), scoring_args={}, suppress_warnings=True)


reportArima(arimaModel, 'Deseasonal', calcPredTestseriesDeseasonal)


%%time
data = seriesBoxCox.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = True)
arimaModel

Wall time: 1min

ARIMA(order=(3, 1, 1), scoring_args={}, seasonal_order=(2, 0, 2, 7),
      suppress_warnings=True, with_intercept=False)


reportArima(arimaModel, 'BoxCox', calcPredTestBoxCox)


%%time
data = seriesDiffHistory.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = False)
arimaModel

C:\Users\herik\AppData\Roaming\Python\Python38\site-packages\pmdarima\arima\_validation.py:62: UserWarning:

m (7) set for non-seasonal fit. Setting to 0

Wall time: 11.6 s

ARIMA(order=(3, 0, 3), scoring_args={}, suppress_warnings=True)


reportArima(arimaModel, 'Stationary', calcPredTestStationary)


result


# Redefinindo variaveis globais para o treino utilizando LSTM

BATCH_SIZE = 30


# hypertuning do algoritmo de LSTM
def GSLSTM(data, learning_rate, calcFunction):
    mse = []
    
    # Crossvalidation para cada parte do conjunto
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]

        try:
            with catch_warnings():
                filterwarnings('ignore')
                
                # Normalização e reshape do conjunto de treino
                ss = StandardScaler()
                ss.fit(train.values.reshape(-1, 1))
                train_input = ss.transform(train.values.reshape(-1, 1))
                
                # Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
                test_input = train_input[-(FEATURES + 1):]
                test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
                train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)
                
                # Reset da sessão
                clear_session()
                
                # Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
                initializer = GlorotUniform(seed = SEED)
                model = Sequential()
                
                # 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
                model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
                model.add(Dense(1, kernel_initializer = initializer))
                model.add(BatchNormalization())
                model.add(Dense(1, kernel_initializer = initializer))
                
                # Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
                early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
                
                # Treinando o modelo com otimizador Adam
                model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
                h = model.fit(train_gen, epochs = EPOCHS, callbacks = [early_stopping], verbose = False)
                pred = np.empty(test.shape[0])
                
                # Realizando predições no conjunto de teste
                for i in range(len(pred)):
                    prediction = model.predict(test_gen, verbose = False)
                    pred[i] = prediction
                    test_input = np.delete(test_input, 0, 0)
                    test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
                    test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
                
                # Reorganizando o shape e chamando a função de calculo
                pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
                pred, test, _ = calcFunction(train, pred, test)      

                mse.append(mean_squared_error(pred, test))
                
        except:
            mse.append(-1)
        
    return np.mean(mse)


# Report do algoritmo LSTM
def reportLSTM(data, modelName, calcFunction, learning_rate):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - Long Short Term Memory (LSTM)'
    
    # Crossvalidation para cada parte do conjunto
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
    
        # Normalização e reshape do conjunto de treino
        ss = StandardScaler()
        ss.fit(train.values.reshape(-1, 1))
        train_input = ss.transform(train.values.reshape(-1, 1))

        # Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
        test_input = train_input[-(FEATURES + 1):]
        test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
        train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)

        # Reset da sessão
        clear_session()
        
        # Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
        initializer = GlorotUniform(seed = SEED)
        model = Sequential()
        
        # 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
        model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
        model.add(Dense(1, kernel_initializer = initializer))
        model.add(BatchNormalization())
        model.add(Dense(1, kernel_initializer = initializer))
        
        # Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
        early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
        
        # Treinando o modelo com otimizador Adam
        model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
        h = model.fit(train_gen, epochs = EPOCHS, callbacks = [early_stopping], verbose = False)
        pred = np.empty(test.shape[0])

        # Realizando predições no conjunto de teste
        for i in range(len(pred)):
            prediction = model.predict(test_gen, verbose = False)
            pred[i] = prediction
            test_input = np.delete(test_input, 0, 0)
            test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
            test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)

        # Reorganizando o shape e chamando a função de calculo
        pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Gerando tabela de hypertunning com taxas de learning_rate
learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam.head()


%%time
if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSLSTM(seriesHistory.copy(), x.learning_rate, calcPredTestOriginal), axis = 1)

Wall time: 3min 51s


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()


reportLSTM(seriesHistory.copy(), 'Original', calcPredTestOriginal, 0.0001)


learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam.head()


%%time

if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSLSTM(seriesDeseasonal.copy(), x.learning_rate, calcPredTestseriesDeseasonal), axis = 1)

Wall time: 4min 1s


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


reportLSTM(seriesDeseasonal.copy(), 'Deseasonal', calcPredTestseriesDeseasonal, 0.01)


learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam.head()


%%time

if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSLSTM(seriesBoxCox.copy(), x.learning_rate, calcPredTestBoxCox), axis = 1)

Wall time: 4min 38s


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


reportLSTM(seriesBoxCox.copy(), 'BoxCox', calcPredTestBoxCox, 0.001)


learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam.head()


%%time

if EXECUTE_GRID_SEARCH:
    hyperparam['mse'] = hyperparam.apply(lambda x: GSLSTM(seriesDiffHistory.copy(), x.learning_rate, calcPredTestStationary), axis = 1)

Wall time: 4min 40s


if EXECUTE_GRID_SEARCH:
    display(hyperparam.query('mse == mse.min() and mse != -1'))


reportLSTM(seriesDiffHistory.copy(), 'Stationary', calcPredTestStationary, 0.00001)


result


# Tratando nomes e criando colunas de controle para os resultados gerados
topResult = (
    result 
    .assign(Full_Name = lambda x: x.Algorithm.apply(lambda a: a.split('(')[0]
                                                   .replace('ARIMA', 'Auto Arima')
                                                   .replace('Long Short Term Memory', 'LSTM')))
    .assign(Data_Category = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[0]))
    .assign(Algorithm = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[1].split('(')[0]
                                                   .replace('ARIMA', 'Auto Arima')
                                                   .replace('Long Short Term Memory', 'LSTM')))
    .sort_values('MSE')
)


topResult.head()


# Plot dos resultados obtidos por ordem ascendente do MSE

colors = {'Time Series Regression':'red',
          'Exponential Smoothing':'orange',
          'Auto Arima': 'green',
          'LSTM ': 'blue'}

# plotly figure
fig = go.Figure(layout = go.Layout(yaxis = {'type': 'category', 'title': 'Algoritmo e Categoria'},
                        xaxis = {'title': 'MSE'},
                        title = 'MSE por Algoritmo e Tipo de Dado'))

for t in topResult['Algorithm'].unique():
    topResultFiltered = topResult[topResult['Algorithm']== t].copy()
    fig.add_traces(go.Bar(x = topResultFiltered['MSE'], y = topResultFiltered['Full_Name'], name = str(t),\
                          marker_color = str(colors[t]), orientation = 'h',
                          text = round(topResultFiltered['MSE']), textposition = 'outside'))
    
    
fig.update_layout(yaxis = {'categoryorder':'total descending'}, autosize = False,
                  width = 1450,
                  height = 800)    
    
fig.show()


# Alocando melhor modelo a memória e separando base de treino
data = seriesHistory.copy()
train = data[datetime.date(2017, 1, 1): datetime.date(2018, 6, 16)]


# Treinando modelo baseado dos parâmetros descobertos na fase de modelagem
ES = (
    ExponentialSmoothing(train, trend = 'add', seasonal = 'add', seasonal_periods = FEATURES)
    .fit(remove_bias = True, smoothing_level = 0.2, smoothing_trend = 0.1, smoothing_seasonal = 0.7, method = 'Powell')
)


# Calculando a previsão até o final do ano de 2018
pred = ES.predict(str(data.index[0]), '2018-12-31')


# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = data.index,
                        y = data,
                        name = 'Real'),
             go.Scatter(x = pred.index,
                        y = pred,
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Deseasonal - Exponential Smoothing')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-10-17', 
              annotation_text = 'Previsão base de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'red', opacity = 0.2, line_width = 0)

fig.add_vrect(x0 = '2018-10-17', x1 = '2018-12-31', 
              annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'green', opacity = 0.2, line_width = 0)

pyoff.iplot(fig)


ES.summary()

	index	order_purchase_timestamp
0	2016-09-15	1
1	2016-09-16	0
2	2016-09-17	0
3	2016-09-18	0
4	2016-09-19	0
...	...	...
709	2018-08-25	69
710	2018-08-26	73
711	2018-08-27	66
712	2018-08-28	39
713	2018-08-29	11

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value
0	Original - Time Series Regression	5296.329792	60.784585	42.960631	179.560484	168.008596
1	Deseasonal - Time Series Regression	5081.777760	58.208746	40.122247	179.560484	164.670085
2	BoxCox - Time Series Regression	5328.602465	61.060147	43.254815	179.560484	165.399796
3	Stationary - Time Series Regression	6057.828087	64.323443	46.538896	179.560484	184.267746

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value
0	Original - Time Series Regression	5296.329792	60.784585	42.960631	179.560484	168.008596
1	Deseasonal - Time Series Regression	5081.777760	58.208746	40.122247	179.560484	164.670085
2	BoxCox - Time Series Regression	5328.602465	61.060147	43.254815	179.560484	165.399796
3	Stationary - Time Series Regression	6057.828087	64.323443	46.538896	179.560484	184.267746
4	Original - Exponential Smoothing	4862.574453	57.625764	38.934719	179.560484	173.334854
5	Deseasonal - Exponential Smoothing	4893.456462	57.825251	39.205131	179.560484	172.837068
6	BoxCox - Exponential Smoothing	50045.745866	159.172302	132.177154	179.560484	205.644407
7	Stationary - Exponential Smoothing	35949.069517	154.499151	101.359143	179.560484	174.147463
8	Original - ARIMA(1,1,2)(0,0,2)[7] intercept	11600.540620	89.231877	70.624151	179.560484	194.679166
9	Deseasonal - ARIMA(2,1,1)(0,0,0)[0] intercept	9591.508344	77.052176	59.641647	179.560484	191.565756
10	BoxCox - ARIMA(3,1,1)(2,0,2)[7]	5548.676257	61.910305	43.954604	179.560484	176.695896
11	Stationary - ARIMA(3,0,3)(0,0,0)[0] intercept	6066.983553	65.860816	47.039531	179.560484	163.822791

	learning_rate
0	0.00001
1	0.00010
2	0.00100
3	0.01000
4	0.10000

	learning_rate
0	0.00001
1	0.00010
2	0.00100
3	0.01000
4	0.10000

1. Preparando Conjunto de Pedidos¶

2. Iniciando Analise Seasonal¶

2.1 Teste Estacionário ADF¶

3. Modelagem¶

3.1 TSR (Time Series Regression)¶

3.2 STL seriesDeseasonal - TSR (Time Series Regression)¶

3.3 BoxCox - TSR (Time Series Regression)¶

3.4 Stationary - TSR (Time Series Regression)¶

3.5 Exponential Smoothing¶

3.6 seriesDeseasonal - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesHistory para seriesDeseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.7 BoxCox - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesHistory para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.8 Seasonal - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesHistory para seriesDiffHistory, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.9 ARIMA¶

3.10 seriesDeseasonal - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada de seriesHistory para seriesDeseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.11 BoxCox - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada de seriesHistory para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.12 Stationary - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada original para stationary, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.13 LSTM¶

3.14 seriesDeseasonal - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para seriesDeseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.15 BoxCox - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.16 Stationary - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para stationary, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

4. Comparação¶

4. Previsões Futuras¶

	trend	season	periods	bias	method
0	add	add	7	True	L-BFGS-B
1	add	add	7	True	TNC
2	add	add	7	True	SLSQP
3	add	add	7	True	Powell
4	add	add	7	True	trust-constr

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value	Full_Name	Data_Category
4	Exponential Smoothing	4862.574453	57.625764	38.934719	179.560484	173.334854	Original - Exponential Smoothing	Original
5	Exponential Smoothing	4893.456462	57.825251	39.205131	179.560484	172.837068	Deseasonal - Exponential Smoothing	Deseasonal
1	Time Series Regression	5081.777760	58.208746	40.122247	179.560484	164.670085	Deseasonal - Time Series Regression	Deseasonal
13	LSTM	5293.714092	58.660983	40.197495	179.560484	165.360380	Deseasonal - LSTM	Deseasonal
0	Time Series Regression	5296.329792	60.784585	42.960631	179.560484	168.008596	Original - Time Series Regression	Original

Dep. Variable:	order_purchase_timestamp	No. Observations:	532
Model:	ExponentialSmoothing	SSE	2219272.189
Optimized:	True	AIC	4456.777
Trend:	Additive	BIC	4503.820
Seasonal:	Additive	AICC	4457.479
Seasonal Periods:	7	Date:	Sat, 01 Oct 2022
Box-Cox:	False	Time:	15:51:26
Box-Cox Coeff.:	None

	coeff	code	optimized
smoothing_level	0.2000000	alpha	False
smoothing_trend	0.1000000	beta	False
smoothing_seasonal	0.7000000	gamma	False
initial_level	0.5877576	l.0	True
initial_trend	1.3908240	b.0	True
initial_seasons.0	-2.1075750	s.0	True
initial_seasons.1	-3.1507588	s.1	True
initial_seasons.2	-4.5117666	s.2	True
initial_seasons.3	-6.4540927	s.3	True
initial_seasons.4	18.982655	s.4	True
initial_seasons.5	-2.3152545	s.5	True
initial_seasons.6	-2.2282936	s.6	True

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4